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A)  GENERAL  CONSIDERATIONS 

The  interaction  between  a  molecule  and  a  solid  surface  can  lead  to  a  great  variety  of  elementary  processes 
such  as  elastic,  inelastic  and  reactive.  Of  particular  importance  is  the  dissociative  chemisorption  of  a 
diatomic  molecule 

02  (v,  j)  +  wall  —  >  [02  (v',  j')\ad  *  wall  —  >0  +  0  +  wall  (1) 

where  a  molecule  chemisorbed  at  the  surface  in  a  specific  roto-vibrational  state  (vj)  dissociates  with  the 
two  atoms  adsorbed  or  scattered  into  the  gas-phase.  Of  great  importance  is  also  the  atom  recombination  on 
surfaces: 

0  +  0  +  wall  — »  02  (v,  j)  +  wall  (2) 

here  two  atoms  recombine  thus  forming  a  diatomic  molecule  that  can  be  either  chemisorbed  or  reflected  in 
the  gas-phase  in  a  given  internal  energy  state.  Reactions  (l)-(2)  are  very  often  the  rate  determining  step  of 
complex  heterogeneous  systems  of  interest  in  different  branches  of  industrial  and  technological 
applications,  as  for  example  in  the  ammonia  synthesis,  hydrocarbon  production,  chemical  vapour 
deposition,  etching  and  thin  solid  film  deposition  via  plasma,  nuclear  rector  technologies  [1]. 

Both  processes  are  of  central  importance  in  aero  thermodynamics  and  the  chemistry  of  interstellar  media. 
Thus,  the  recombination  of  atomic  O  and  N  on  silica  and  UHTC  materials  plays  a  central  role  for  the 
thermal  protection  system  of  the  space  shuttles  entering  into  the  terrestrial  atmosphere,  whereas  the 
recombination  of  hydrogen  atoms  on  ice  grains  covered  by  carbon  is,  very  likely,  the  main  source  of 
molecular  hydrogen  observed  in  the  interstellar  media. 

The  interaction  of  chemical  species  with  surfaces  can  lead  to  other  non-reactive  chemico-physical 
processes  such  as  the  inelastic  processes  and  adsorption.  The  adsorption  processes  occur  when  the  particle 
is  trapped  in  a  chemisorption  site  and  its  available  energy  is  not  enough  to  escape  from  the  chemisorption 
potential  well.  The  inelastic  processes  can  be  of  two  types:  direct  and  indirect.  In  the  first  case  the 
molecule  hits  the  surface  and  it  is  scattered  in  the  gas-phase  after  few  bounces.  Due  to  the  interaction,  the 
internal  energy  state  of  the  molecule  leaving  the  surface  is  different  with  respect  to  its  state  before  the 
collision. 

In  the  indirect  scattering,  or  adsorption/desorption  collisions,  the  interaction  occurs  on  a  longer  time  scale 
and  the  collision  is  more  involved:  in  fact,  the  molecule  hits  the  surface  several  times,  it  is  temporarily 
adsorbed  eventually  forming  a  surface  activated  complex  before  desorbing  in  the  gas-phase.  The  energy 
distribution  of  the  scattered  molecules  is  different  according  to  the  different  mechanism  followed  in  the 
interaction. 
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Molecular  Dynamics  Simulations  of  Surface  Processes: 

Oxygen  Recombination  on  Silica  Surfaces  at  High  Temperature 


A  further  class  of  surface  processes  are:  surface  oxidation,  chemical  and  physical  sputtering,  ion 
implantation,  etc.  In  these  cases  atomic  and  molecular  radicals  react  with  the  atoms  of  the  substrate,  thus 
forming  volatile  compounds.  These  processes,  not  discussed  in  this  lecture,  are  effective  in  the  energy 
regime  of  several  eV,  or  KeV,  that  is  higher  compared  to  the  collisional  energy  regime  of  the  processes 
(l)-(2). 

The  dynamical  quantity  typically  associated  to  the  process  (1)  is  the  state-to- state  sticking  coefficient  S0 
which  corresponds  to  the  dissociation  probability.  So  is  a  function  of  the  several  variables:  So=  So 
(Ekin\{9i,  (pi}\v,j,Ts ),  where  Ekin  is  the  impact  kinetic  energy,  v  and  j  the  vibrational  and  rotational  quantum 
number  of  the  incident  molecule  respectively,  Ts  is  the  surface  temperature.  { 9t,  cpj  are  the  angles  that 
define  the  relative  orientation  of  the  molecule  with  respect  to  the  surface  plane  (  the  collisional  system, 
sketched  in  Fig.l,  is  described  in  a  cartesian  frame  of  reference  having  the  Z  axis  orthogonal  to  the 
surface  plane  and  the  (X,Y)  axis  laying  on  the  surface  top  layer). 


Fig.  1 :  The  collisional  system  is  sketched.  The  dynamics  is  described  with  respect 
to  a  cartesian  frame  of  reference  having  the  Z  axis  orthogonal  to  the  surface 
plane  and  the  (X,Y)  axis  laying  on  the  surface  top  layer. 

Different  averaged  quantities  can  be  obtained  from  the  state-to- state  sticking  coefficient.  Thus  the  state- 
average  sticking  coefficient  <So(Ekin\{9iy  (pJ\Ts)>,  usually  measured  in  molecular  beam  experiments,  is 
obtained  by  averaging  So  over  a  known,  not  necessarily  Boltzmann  population  distribution  for  the  v  and  j 
states.  The  statistical  rate  constant  K(TS)  used  in  kinetic  modelling  can  be  obtained  by  further  averaging  So 
over  the  incident  angles  and  the  Maxwell  energy  distribution  function  of  the  gas-phase  molecule.Then  we 
get  K(Ts)=Kcoll*<So(Ts)>. 

In  addition  to  S0 ,  other  collisional  data  needed  to  characterize  the  surface  reactions  and  to  include  them  in 
kinetic  modelling  are:  the  energy  accommodation  coefficient  (3(Ts),  the  roto-vibrational  distribution  of  the 
scattered  particles  N(v,jlTs),  the  angular  distribution  in  intensity  and  energy  I(0,cplTs),  the  lifetimes  of  the 
adsorbed  species. 

Collisional  data  can  be  determined  by  using  a  large  variety  of  experimental  techniques  that  have  been 
developed  to  probe  directly  or  indirectly  surface  processes.  Nevertheless,  the  observation  of  surface 
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processes  is  not  an  easy  task  and  in  fact  it  is  a  matter  of  fact  that  collisional  data  for  heterogeneous 
processes  are  very  sparse  (if  not  absent)  despite  the  large  amount  of  information  emerged  from  different 
experimental  works.  Rate  coefficients  are  still  poorly  known  and  somehow  dubious,  particularly  for  wall 
processes  occurring  under  extreme  thermal  conditions  as  those  met  in  aerothermodynamics. 

In  Fig.  2  some  of  the  experimental  results  reported  in  the  literature  for  the  recombination  coefficient  y  for  O 
atoms  recombination  on  silica-based  materials  are  reported  in  a  large  range  of  the  surface  temperature  [2-4]. 
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Fig.  2:  The  recombination  coefficient  y  for  O  atom  recombination  on  silica-based 
materials  is  reported  as  a  function  of  the  surface  temperature. 


In  the  high  temperature  regime  the  results  are  somehow  contrasting,  also  due  to  the  large  experimental 
uncertainties.  From  the  theoretical  point  of  view  various  kinetics  schemes  have  been  developed  to  describe 
the  kinetics  of  this  reaction  [5,6,7]. 

In  these  studies  some  physical  quantities,  typically  the  activation  energy,  surface  coverage,  the  number  of 
active  sites  and  the  so-called  ’steric  factor',  are  treated  as  unknown  parameters  and  varied  to  fit  the 
experimental  data. 

Due  to  the  lack  of  accurate  collisional  data,  a  complete  understanding  of  elementary  surface  processes 
involving  O,  N,  C,  CO  on  silica-based  materials  has  not  yet  completely  achieved  and  a  number  of 
fundamental  aspects  are  open  to  question.  Among  the  others,  the  most  critical  aspects  discussed  in 
molecule-surface  interaction  studies  concern,  on  the  one  hand  the  internal  state  selectivity  of  the 
collisional  coefficient  and,  on  the  other  hand  the  energy  transfer  mechanisms  that  control  the  surface 
processes.  Further  elements  of  complexity  are  introduced  by  the  surface  structure  and  surface  coverage 
effects  on  the  catalytic  activity  of  the  substrate. 

Possible  answers  to  these  issues  can  be  given  by  performing  Molecular  Dynamics  (MD)  calculations 
[8,9,10].  The  MD  approach  has  some  important  features  in  that  predicts  quantities  that  cannot  be  easily 
measured  in  experiments:  quite  often,  the  state- selected  sticking  coefficient  So,  or  the  energy 
accommodation  coefficient  p  can  be  established  only  on  a  theoretical  ground. 
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The  MD  simulation  of  surface  chemical  processes  is  quite  complex  both  from  the  formal  and  the 
computational  point  of  view.  The  multidimensional  nature  of  the  collisional  system,  and  thus  the  large 
number  of  degrees  of  freedom  introduced  by  the  presence  of  the  solid,  increases  considerably  the 
complexity  of  the  dynamics.  As  a  consequence,  in  the  last  few  years  a  variety  of  collisional  methods  [8,9] 
have  been  developed  to  simulate  reactive  surface  processes,  ranging  from  ab  initio  quantum  mechanics 
calculations  of  reduced  dimension,  semiclassical  and  mixed  semiclassical-wave  packet  propagation 
techniques,  statistical  simulations,  classical  and  quasi-classical  methods.  MD  simulations  imply  the 
solution  of  three  main  problems:  i.)  the  first,  probably  the  most  critical  one  that  would  prevent  an  accurate 
description  of  the  surface  processes  concerns  the  determination  of  the  potential  energy  surface  (PES) 
where  the  dynamics  takes  place,  that  is  the  determination  of  the  forces  exerted  between  the  atoms  in  the 
gas-phase  and  the  atoms  of  the  surface,  ii.)  The  second  concerns  the  building  up  of  a  ‘model’  sample  of 
the  solid  substrate,  iii.)  Finally  the  third  aspect  to  solve  concerns  the  development  of  a  suitable  collisional 
model  to  follow  the  nuclear  motions  of  the  particles  propagating  in  the  gas  phase  and  interacting  with  the 
solid  surface. 

In  the  following  sections  some  interesting  results  obtained  from  MD  simulations  of  atom  recombination 
processes  over  model  silica  surfaces  that  could  be  of  interest  in  aero  thermodynamics  are  briefly  reviewed 
and  discussed. 

B)  INTERACTION  POTENTIAL  DETERMINATION  FOR  0,N/SILICA 

B.l)  Ab  initio  methods 

Determining  the  PES  where  the  reaction  takes  place  is  a  prerequisite  for  any  MD  calculations.  In  principle, 
this  issue  is  addressed  using  computational  methods  developed  in  quantum  chemistry  [10].  Among  the 
proposed  approximate  schemes,  the  most  promising  approach  is  the  density  functional  theory  (DFT) 
[9,11-13],  which  is  accurate  near  the  chemisorption  well,  but  less  accurate  for  long-range/weakly  bound 
interactions.  The  accurate  determination  of  the  interaction  potential  for  heterogeneous  systems  poses 
different  problems,  both  from  the  theoretical  and  the  computational  point  of  view.  Assuming  that  the 
relaxation  of  the  lattice  atoms  can  be  neglected,  the  PES  is,  in  principle,  a  multidimensional  function 
which  depends  for  diatom- surface  interaction,  upon  six  coordinates,  that  is  the  surface  is  too  much 
complex  and  the  number  of  interactions  is  too  large  for  ab  initio  electronic  structure  calculations,  either 
cluster  or  slab  calculations,  involving  a  large  number  of  atoms  and  electrons.  As  a  consequence,  also  due 
to  the  large  computational  time  demanded  to  assure  the  convergence  of  calculated  energies,  complete  fully 
dimension  electronic  structure  calculations  have  been  performed  for  few  elementary  systems,  notably 
those  involving  H2,  CO,  02  and  others  diatomics  on  noble  and  transition  metals.  The  'true’  interaction 
potential  remains  rather  unknown  for  a  large  class  of  heterogeneous  systems,  in  particular  for  the  catalytic 
systems  of  interest  in  aerothermodynamics.  In  fact  uncertainties  exist  on:  the  adsorption  energies  and  their 
dependence  on  the  chemisorption  site,  the  size  and  the  nature  of  the  activation  energies,  the  energy 
barriers  to  surface  diffusion. 

Due  to  the  lack  of  ab  initio  calculations  one  has,  quite  often,  to  rely  on  low-dimension  ab  initio  ‘model’ 
potentials  where  the  dependence  of  the  interaction  potential  upon  the  intramolecular  distance  of  the  diatom 
and  the  distance  from  the  surface  is  considered.  PES  of  reduced  dimension  can  be  useful  for  understanding 
the  qualitative  aspects  of  dissociation/recombination  processes,  but  for  a  quantitative  description  a  PES 
that  includes  the  dependence  on  the  two  (X,Y)  surface  plane  coordinates  is  necessary. 

An  accurate  two-dimensional  adiabatic  PES  for  N,N2  interacting  with  a  silica  surface  has  been  recently 
calculated  [14]  using  the  size-scalable  cluster  approach  with  SixOy  clusters  of  increasing  size  cleaved  from 
the  (3-cristobalite  unit  cell.  In  this  study  the  hybrid  Hartree-Fock  Self  Consistent  Field  (HF-SCF)  and  the 
DFT  method  was  applied  and  the  B3LYP  functions  [15]  used  in  order  to  take  into  account  the  exchange 
and  the  electron  correlation  contribution  in  the  interaction.  The  interaction  potential  is  calculated  by 
keeping  the  Si02  geometry  fixed  at  the  experimental  values  of  the  P-cristobalite  cell  and  extending  the 
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N/N2-SixOyHz  cluster  model  in  the  intervals  (x=l,3,7;  y=2,4,6,14;  z=2,6,14).  All  the  calculations  were 
carried  out  with  the  Gaussian  03  package  [16].  Complete  scans  of  the  RN.Si  distance  in  the  whole  range  of 
energy  interaction  (Rn_Si  =  0.90  A,  RN_Si=4.00  A)  for  just  the  perpendicular  configuration  (that  is,  with  the 
N-surface  angles  fixed  to  90°)  were  performed.  The  N-N  distance  used  in  all  N2-surface  PES  calculations 
was  kept  fixed  to  1.10  A,  close  to  the  spectroscopic  value  [17].  As  far  as  the  N-SixOyHz  interaction  is 
concerned,  both  doublet  and  quartet  total  electronic  spin  states  were  considered  along  the  PES  scan,  since 
the  electronic  spin  state  induced  by  chemisorption  is  not  known.  The  calculated  binding  energies  are 
reported  in  Fig.  3.  As  expected,  N  is  chemisorbed  with  a  binding  energy  of  about  2.79eV  at  a  minimum  N- 
Si  distance  around  1.65  A.  N2  is  slightly  physisorbed. 


Fig.  3:  Full  line:  the  interaction  potential  for  N  interacting  with  S^C^He  cleaved  from  the 
p-cristobalite  surface  is  reported  as  a  function  of  the  distance  of  the  N  atom  from  the 
Si  active  atom.  The  dashed  line  shows  the  interaction  potential  of  N2(  req) 
approaching  the  silica  cluster  in  the  perpendicular  configuration  on  top 
of  the  Si  active  atom.  Results  obtained  at  the  DFT-B3LYP/6-311+G*  level. 


From  the  calculated  quantum  mechanical  interaction  energies  a  LEPS-type  interaction  potential  can  be 
obtained  [18,19].  In  order  to  account  for  the  motion  of  the  lattice  atoms  explicitly  considered  in  the 
collision  dynamics,  the  LEPS  potential  is,  in  this  case,  obtained  as  a  sum  of  pair-wise  N/silica  atom 
interaction  and  it  is  given  by: 


U,(r)  +  U2(R.,  )  +  Us(Ru  )-Jtt(r)  +  (A,(R'J)  +  A,(Ru  )?  -A,(rfa(R„  )+A,(Su  >)  (3) 


where  the  sum  is  over  the  lattice  atoms,  W  being  the  total  number  of  atoms  in  the  silica  surface,  r  is  the  N- 
N  interatomic  distance,  Ra  i  and  Rb  i  are  the  distances  of  nitrogen  atom  a  and  nitrogen  atom  b  from  the  i-th 
atom  in  the  lattice,  respectively.  Uk  ( k=l,3 )  and  Ak  ( k=l,3 )  are  related  to  the  two-body  Coulomb  and 
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exchange  integrals,  respectively,  in  the  energy  expression  for  the  interaction  of  the  solid  surface  with  the 
molecule  [18,19]  and  are  given  as  a  combination  of  modified  Morse  functions  with  the  two  additional  Sato 
parameters.  The  corresponding  potential  parameters  are  obtained  by  fitting  the  full  potential  given  by  Eq. 
(3)  to  the  ab  initio  DFT-B3LYP  energies  (Fig. 3).  Since  the  determination  of  the  activation  barrier  for  the 
N2  formation  would  imply  massive  multidimensional  DFT  calculations,  planned  for  the  near  future,  the 
Sato  parameters  are  taken  as  a  free  parameters  and  changed  in  order  to  get  PES’s  with  different  barriers  in 
the  recombination  reaction  channel.  This  allows  to  study  the  sensitivity  of  the  calculated  collisional 
reaction  coefficients  upon  the  reaction  barrier  height. 

In  Fig.  4  the  potential  contour  map  of  the  PES  for  N2  interacting  in  the  perpendicular  geometry  with  the 
active  Si  surface  atom  is  shown.  The  interaction  potential  is  plotted  as  a  function  of  the  distance  from  the 
N  atom  closest  to  the  surface,  and  the  N-N  bond  distance.  The  PES  shown  is  obtained  assuming  a  value  of 
the  Sato  parameters  that  gives  zero  activation  energy. 


Fig.  4:  PES  for  N2  interacting  in  the  perpendicular  geometry  with  the  active  Si  surface  atom. 
The  interaction  potential  (in  eV)  is  plotted  as  a  function  of  the  distance  from 
the  N  atom  closest  to  the  surface  and  the  N-N  bond  distance. 


The  DFT  cluster  approach  has  been  recently  extended  to  the  0,02/ZrB2  system  and  preliminary  results  for 
O  chemisorbed  on  top  of  a  Zr  atom  have  been  obtained  by  us. 

B.2)  Approximate  semiempirical  methods 

For  molecule-surface  systems  potential  surfaces  can  be  inferred  following  a  semiempirical  approach  by 
using  the  spectroscopic  data  available  for  the  chemisorbed  species,  i.e.  equilibrium  geometries,  adsorption 
energies,  inter  and  intra-molecular  vibrational  frequencies.  The  idea  is  to  start  with  a  parametric 
interaction  potential  and  to  calculate  the  potential  parameters  such  that  the  potential  is  able  to  reproduce 
all  the  observed  spectroscopic  constants.  This  approach  has  been  followed  to  derive  a  realistic  interaction 
potential  for  CO/Pt(l  1 1)  [20],  a  system  considered  as  a  prototype  in  heterogeneous  catalysis. 
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Unfortunately,  due  to  the  lack  of  spectroscopic  data  for  the  0,02/silica  system,  this  approach  cannot  be 
followed.  However  a  tentative  interaction  potential  can  be  constructed  semiempirically  on  the  ground  of 
reliable  physical  considerations.  The  method,  developed  to  study  the  interaction  of  O  with  two  silica 
polymorphs  [21,22],  (3-cristobalite  and  (3-quartz,  although  quite  crude  has  been  successful  in  predicting  the 
recombination  coefficient  y  for  the  O+O/silica  system  at  high  temperature.  As  usual,  the  O-silica  potential 
is  expressed  as  a  sum  of  pair-wise  interactions,  assuming  that  the  Ogas-Si  interaction  is  the  same  as  the  O- 
Si  interaction  in  the  lattice,  while  the  interaction  between  the  O  atom  in  the  gas-phase  and  the  O  of  the 
quartz  surface  is  assumed  as  that  known  for  the  0-0  interaction  in  the  gas-phase.  The  validity  of  this  latter 
assumption  is  discussed  later  in  the  paper  .Thus  we  have: 

V0.s=Z  SV1KS(R,)t|vM(R,.) 

a=l  ,2  _  j=l  j  =1 

where  the  first  term  is  the  BKS  potential  [23]  for  the  interaction  of  the  gaseous  O  atoms  and  the  Si  lattice 
atoms.  Moreover,  since  the  BKS  potential  at  O-Si  distances  shorter  than  1.5  A  is  attractive,  in  order  to  get 
the  expected  correct  behaviour,  the  following  potential  is  taken  for  the  Ogas-Si  interaction: 


+  Vr 


(4) 
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O-Si 
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R O-Si  ^0 
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(5) 


where  the  repulsive  and  the  attractive  part  of  the  BKS  have  been  retained  (and  the  Coulomb  charge-charge 
interaction  dropped  out).  An  appropriate  switching  function  has  been  introduced  in  order  to  smooth  the 
repulsive  potential  to  the  C6  attractive  interaction  at  the  larger  distances.  R0  was  set  to  2.5  A.  The  second 
term  in  Eq.  (4)  is  the  interaction  potential  between  the  gaseous  O  atoms  and  the  O  lattice  atoms  and  it  is 
given  as: 


Vq-o  ~  A0~o  e 


(6) 


where  the  parameters  A 0-o  and  b0_ 0  are  those  known  from  the  0-0  interaction  in  the  gas-phase  [21]. 
Finally,  the  last  term  in  Eq.  (4)  represents  the  Coulomb  interaction  between  the  charge  on  the  Si  and  O 
lattice  atoms  and  the  charges  of  the  molecular  oxygen  due  to  the  quadrupole  charge  separation  in  02. 
Further  details  about  these  two  last  terms  of  the  potential  can  be  found  in  Ref.  [21]. 


In  Fig. 5,  the  two-dimensional  interaction  potential  is  shown  for  02  interacting  in  the  perpendicular 
orientation  on  top  of  a  Si  atom.  The  potential  is  plotted  as  a  function  of  the  0-0  internuclear  distance  and 
the  distance  of  the  closest  O  to  the  surface.  In  this  configuration,  the  PES  shows  that  the  02  formation 
occurs  without  activation  energy  so  that  this  configuration  is  very  likely  the  most  effective  for  the 
recombination  process. 
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Fig.  5:  Two-dimensional  interaction  potential  for  O2  interacting  in  the  perpendicular  orientation 
on  top  of  a  Si  atom  of  a  p-quartz  surface.  The  potential  (in  eV)  is  plotted  as  a  function  of 
the  0-0  internuclear  distance  and  the  distance  of  the  closest  O  atom  to  the  surface. 


C)  RECOMBINATION  DYNAMICS  OF  O  ON  ‘MODEL’  SILICA  SURFACES 
C.l)  The  semiclassical  collisional  method 

In  this  section  two  aspects  relevant  to  the  MD  approach  to  catalysis  are  discussed:  firstly  we  briefly 
discuss  the  semiclassical  collisional  model  developed  to  describe  the  dynamics  of  the  O  atom 
recombination  and  the  02  dissociation  on  silica,  secondly  the  static  properties  of  a  silica  surface  are 
examined.  Compared  to  a  fully  classical  MC  simulation,  a  much  detailed  and  physically  coherent 
description  of  the  surface  phenomena  can  be  reached  using  the  semiclassical  time-dependent  collisional 
method  [24].  This  method  is  able  to  describe  the  most  fundamental  features  of  the  molecule/surface 
interaction:  desorption  (inelastic  scattering)  and  sticking,  diffusion,  dissociation  and  atom  recombination. 
The  following  aspects  can  be  incorporated  in  the  model:  lattice  geometry,  force  constants,  electronic 
structure  (band  gap,  Fermi  level,  etc). 

The  method  is  semiclassical  in  that,  while  the  nuclear  motions  of  the  atoms/molecules  approaching  the 
surface  are  treated  classically,  the  vibrational  motions  of  the  lattice  atoms  and  the  electrons  are  treated 
quantum-mechanically.  Thus,  the  classical  Hamilton’s  equations  of  motion  for  the  gas-phase  species  and 
the  time-dependent  Schrodinger  equations  of  motion  for  the  phonons  and  the  electrons  of  the  substrate  are 
solved  self-consistently.  The  time  evolution  of  the  phonons/electrons  is  dynamically  coupled  to  the  time 
evolution  of  the  chemical  species  through  the  definition  and  calculation  of  the  effective  Hamiltonian 


(7) 


H  ^  is  contribution  due  to  the  excitation/deexcitation  processes  of  the  surface  phonons  and  it  is  defined 
as  the  expectation  value  of  the  interaction  potential  over  the  total  wave  function  of  the  phonon  state  at  a 
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given  temperature.  Similarly,  the  effective  Hamiltonian  due  to  the  excitation/de-excitation  of  the  electron 


e-h 

eff 


holes  in  the  substrate,  H 
electrons: 

n$ 


,  is  given  by  averaging  the  coulomb  potential  over  the  wave  function  of  the 


=('?«.,  (<KM) 


(8) 


where  t  is  the  interaction  time,  Vml( r,  R)  is  the  interaction  potential  between  the  gas-phase  species  and  the 


surface  atoms  determined  as  described  in  the  previous  section,  |  Wph )  is  the  total  wave  function  of  the 


phonon  state  at  a  given  surface  temperature  Ts ,  Vc(r,R)  is  the  coulomb  potential  due  to  the  interaction 


between  the  electron  charge  in  the  solid  and  the  charge  temporally  transferred  to  the  molecule,  |  ^ e_h )  is 
the  total  wave  function  of  the  electrons  perturbed  by  the  external  coulomb  force. 

The  total  hamiltonian  for  the  molecular  motion  of  a  diatomic  molecule  is  obtained  by  adding  to  the  usual 
hamiltonian  for  the  free  molecule  the  effective  hamiltonian  arising  from  the  molecule- surface  interactions: 

>2 
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i  pA 

=  -Y  ^L  +  V(r)  +  H 


+  H 


e-h 

eff 


+  ^Eph  +  EEe_h 


(9) 


where  P  ia  is  the  a-th  cartesian  component  of  the  momentum  of  atom  i  and  mi  the  corresponding  mass. 
A Eph  and  A Ee_h  are  the  energy  exchanged  with  the  phonons  and  the  electrons  respectively,  V(r)  is  the 
intra-molecular  potential  of  the  gas-phase  molecule. 

The  effective  Hamiltonian  H ^  due  to  the  phonon  interaction  can  be  obtained  within  the  well  known 

linearly  and  quadratically  forced  harmonic  oscillator  model  [24].  This  leads  to  an  analytical  solution  for 
the  time  evolution  of  the  phonon  wave  function  and,  consequently,  for  the  excitation/de-excitation 

amplitudes  of  the  phonons  in  each  normal  modes.  Finally,  an  analytical  expression  of  is  obtained  in 
terms  of  the  ‘phonon  excitation  strengths’  r|k(t): 

H%='Zvl1>(R(t))T7*(t)  (H>) 

k 

where  V(k1}  is  the  linear  external  force  exerted  on  the  £-th  normal  mode,  R(t)  is  the  position  vector  of  the 

incoming  particle  at  t  time.  R(t)  is  the  classical  trajectory  that  is  solution  of  the  Hamilton’s  equations  of 
motion.  r/k(t)  are  given  in  terms  of  the  Fourier  components  of  the  forces  exerted  between  the  gaseous 
species  and  the  surface. 

r  1  d(  A E  h  )  r  1 

r/k(t)  =  -  J  dt'——— - [I  k  cos( (Okt)  +  Isk  sin( cokt)\  (11) 

J  ticok  dpk 

/•+00 

Ic,k  =  dtV(kx>(  R(  t  ))cos(  a>kt )  (12) 

J— oo 

(Ok  are  the  frequencies  of  the  phonon  normal  modes. 

The  method  outlined  here  has  two  important  features:  firstly,  since  the  quantum  phonon  dynamics  is 
solved  exactly  both  the  single  and  the  multi-phonon  processes  are  accounted  for  in  the  collisional  model, 
and  secondly  the  effective  Hamiltonian  depends  explicitly  upon  the  temperature  of  the  solid.  This  later 
effect  is  due  to  the  energy  symmetrization  and  to  the  use  of  a  thermal  phonon  distribution  instead  of  a  pure 
phonon  eigenstate  for  the  initial  phonon  state. 
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A  further  important  feature  is  the  energy  exchanged  between  the  solid  substrate  (phonons  and/or  electrons) 
and  the  chemical  particles.  In  the  case  of  surface  processes  assisted  by  the  phonon  excitation  mechanism, 
the  exchanged  energy  is  given  as: 

k  nk  k 

where  (  En  —E0)  is  the  energy  exchanged  in  the  transition  nk  <—  nk  between  the  quantum 

state  nk  and  nk  of  the  k-th  phonon  mode,  p  0  is  the  Boltzmann  distribution  of  the  phonon  energies, 

V  V  nk 

P  o  is  the  transition  probability. 

nk  ~^nk 

A Ek  is  the  energy  loss  (+)  or  gained  (-)  from  the  molecule  due  to  the  excitation/de-excitation  phonon 
processes  in  the  substrate. 

The  method  is  formally  rather  involved,  but  nevertheless  it  is  computationally  affordable  still  providing  a 
realistic  description  of  the  gas-surface  interactions  in  their  full  dimension.  The  method  has  been  applied  to 
several  heterogeneous  systems,  including  the  CO  and  C02  formation  on  a  platinum  surface  [25,26]  and  the 
inelastic  scattering  of  Xe  atoms  from  a  corrugated  GaSe(OOl)  surface  [27]. 

C.2)  Model  silica  surfaces  and  phonon  dynamics 

In  order  to  perform  MD  calculations  a  sample  of  the  solid  substrate  where  the  surface  reactions  occur  must 
be  built  up  and  incorporated  in  the  collisional  model.  Model  surfaces  can  be  constructed  starting  from  the 
known  arrangement  of  the  Si  and  O  atoms  in  the  unitary  cell  of  a  specific  silica  polymorph  [28],  then  the 
unitary  cell  is  repeated  several  times  in  the  3D  space  in  order  to  have  a  sufficiently  large  sample.  The 
obtained  surface  can  be  considered  as  an  ideal,  perfect  surface.  In  fact,  while  the  bulk  structure  of  real 
surface  is  known,  the  atom  arrangement  on  the  top  most  layers  is  not  due  to  the  covalent  character  of 
silica.  Indeed,  MD  simulations  of  the  static  properties  of  silica  have  revealed  a  very  complex  structure  of 
the  surface  layers  which  appear  reconstructed  with  structural  and  bond  defects(see  Fig.  6)  [29]. 


Two-membered 

rings 


P-cristobalite: 

Si-O  bonding  angle  =  109.48° 
Si-O  bonding  length  =  1.54  A 
0-0  distance  =  2.5  A 


Fig.  6:  A  reconstructed  surface  is  depicted. 
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In  addition  to  that,  the  chemical  and  structural  properties  of  silica  targets  used  in  cell  experiments,  or  in 
real  systems,  can  be  strongly  connected  to  the  chemistry  and  the  physics  of  the  plasma  environment.  As  a 
consequence,  the  surface  reactivity  can  be  different  under  different  operative  plasma  conditions.  Since  the 
real  structure  of  silica  surfaces  cannot  be  easily  modelled,  also  due  to  the  lack  of  structural  data 
determined  from  in  situ  observation,  practically  all  the  MD  simulations  are  carried  out  on  flat  and  perfect 
surfaces.  This  poses  some  problems  when  the  collisional  data  obtained  in  MD  studies  are  compared  with 
those  determined  in  cell  experiments.  On  the  other  hand,  MD  data  can  be  directly  compared  (  and 
validated  )  with  data  measured  in  molecular  beam  experiments,  where  a  molecule/atom  in  a  well  defined 
internal  state  hits  a  well  defined  single  crystal  surface  (as  in  the  MD  simulations). 

In  Fig. 7  the  top  view  of  the  (3-quartz  and  (3-cristobalite  surfaces  used  in  recent  MD  calculations  are  shown. 


Fig.  7:  Top  view  of  the  p-cristobalite  (left)  and  p-quartz  (right) 
surfaces  recently  used  in  MD  calculations. 


In  order  to  describe  the  phonons  dynamics  of  the  silica  surface  a  reliable  interatomic  potential  is 
necessary.  Several  interatomic  potentials  for  silica  polymorphs  are  available  in  the  literature  each  of  them 
being  able  to  reproduce  on  a  quite  good  level  different  structural  properties  of  the  materials.  The  semi- 
empirical  BKS  potential  [23]  is  among  the  prime  potentials  used  for  silica  polymorphs.  The  BKS  potential 
is  a  pair-wise  potential  given  by: 

v,ks(r„  >= v'v'"  ~%+I4A^ir  <14) 

Rij  Rii 

where  the  first  term  describes  the  purely  repulsive  interaction  between  two  lattice  atoms,  i  and  j,  Ry  being 
the  distance  between  the  two  atoms.  The  second  term  is  the  C6  attractive  dispersion  term,  while  the  last 
term  is  the  coulomb  interaction  between  the  charged  Si  and  O  atoms.  Being  the  BKS  potential  a  two-body 
interaction  potential  with  few  parameters,  the  computational  time  needed  to  describe  the  dynamics  of  the 
lattice  phonons  is  lower  than  that  required  by  other  interaction  potentials  and  its  results  in  reproducing 
various  structural  material  properties  are  often  better. 

The  phonon  frequency  spectrum,  {ook},  for  the  assumed  silica  surface  can  be  deduced  from  the  normal 
mode  analysis  of  the  3D  lattice  atom  vibrations.  The  dynamical  matrix  H  of  the  potential  force  constants  is 
constructed  by  calculating  the  second  derivatives  of  the  pair-wise  BKS  potential  applied  to  the  first  nearest 
atoms  of  the  lattice 
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_  Wy 

i<j<N 

dRiadRip 

where  Kiajp  stands  for  the  inter-atomic  force  constants  between  the  two  nearest  lattice  atoms  i  and  j,  a 
and  (3  are  the  components  of  their  respective  cartesian  coordinates.  N  is  the  total  number  of  atoms  in  the 
lattice.  Then  the  secular  equation  is  solved: 

\H-E<jo\  =  Q  (16) 

and  the  eigenvalues  {cok}  corresponding  to  the  frequencies  for  the  localized  vibrational  normal  modes  of 
the  lattice  are  obtained  (E  is  the  unitary  matrix). 

In  Fig.  8,  the  frequency  distribution  calculated  for  a  model  (3-quartz  surface  is  shown  for  a  3D  crystal 
consisting  of  198  atoms  displayed  over  8  layers  [22]. 


Fig.  8:  Frequency  distribution  calculated  for  a  3D  model  crystal 
of  p-quartz  consisting  of  198  atoms  displayed  over  8  layers. 


C.3)  E-R  and  L-H  Oxygen  recombination  dynamics  on  silica 

The  oxygen  atom  recombination  on  surfaces  is  generally  described  according  to  two  extreme  collisional 
mechanisms,  the  so-called  E-R  and  L-H  mechanism.  The  E-R  reaction  is  a  two  steps  process,  firstly  the 
oxygen  is  adsorbed  at  the  surface  in  a  specific  surface  site: 

O  +  silica - >  Oad  +  silica  (17) 

Then  the  adsorbed  oxygen  reacts  with  an  oxygen  atom  approaching  the  surface  from  the  gas-phase. 

O  +  Oad  +  silica - >  02(v,j)  +  silica  (18) 

The  first  step  is  a  thermal  equilibrium  process,  the  second  reaction  is  the  rate-determining  step. 
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According  to  the  L-H  mechanism  the  reaction  occurs  between  two  adsorbed  atoms  which  propagate  on  the 
surface  from  site  to  site  and  then,  eventually,  recombine. 

Oad  +  Oad  +  silica - >  02(v,j)  +  silica  (19) 

MD  simulations  have  been  performed  for  both  reactions  catalysed  by  P-cristobalite  in  the  temperature 
range  from  Ts=600K  to  Ts=1600K,  using  a  semiempirical  interaction  potential  determined  according  to 
the  basic  considerations  discussed  in  a  previous  section.  Details  can  be  found  in  Ref.s  [21,30].  In  Ref.[21] 
the  effect  of  first  layer  composition  on  the  recombination  reaction  was  also  explored.  In  Fig.  9. a  the  E-R 
recombination  probability  at  Ts=1000K  is  reported  for  an  oxygen  covered  and  a  clean  P-cristobalite 
surface,  full  dot  and  open  dot  respectively. 
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Fig.  9:  (a)  The  E-R  oxygen  recombination  probability  at  Ts=1000K  is  reported  for  a  clean  (open  dots) 
and  oxygen  covered  (full  dots)  model  surfaces  of  p-cristobalite.  (b)  The  E-R  recombination 
probability  of  O  on  p-cristobalite  (full  dot)  and  p-quartz  (open  dot)  is  reported  as  a  function 
of  the  impact  energy  of  the  gas-phase  oxygen  atom.  The  surface  temperature  is  Ts=1000K. 


In  Fig. 9b  the  E-R  recombination  probability  of  O  on  P-cristobalite  (full  dot)  and  on  P-quartz  (open  dot)  is 
shown  as  a  function  of  the  impact  energy  of  the  gas-phase  oxygen.  The  surface  temperature  is  Ts=1000K. 

In  the  trajectory  simulations  the  oxygen  atom  sitting  on  the  surface  is  in  thermal  equilibrium  with  the 
surface  and  randomly  adsorbed  within  the  unit  cell  area,  the  initial  Oad-surface  distance  being  2.A.  The 
incoming  gas-phase  O  atom  approaches  the  silica  surface  in  the  perpendicular  orientation,  the  initial  (X,Y) 
position  coordinates  are  randomized  within  the  unit  silica  cell  and  Z  far  in  the  asymptotic  region  (  Z=10  A 
at  t=0 ). 

From  the  calculated  recombination  probabilities  the  recombination  coefficient  y(Ts)  is  calculated 
assuming  that  a  flux  of  oxygen  atoms  hits  the  surface  with  a  maxwellian  kinetic  energy  distribution 
corresponding  to  the  surface  temperature. 

The  y(Ts)  coefficient  is  obtained  by  averaging  the  calculated  recombination  probabilities  over  a  flux  of  O 
atoms  having  a  Maxwellian  kinetic  energy  distribution  at  a  given  temperature  (usually  the  temperature  of 
the  surface): 


Ze-r 
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(Pe-R  (E kin  1 
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kin 


(20) 


where  KB  is  the  Boltzmann  constant. 
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The  L-H  reaction  mechanism  was  also  simulated  and  the  recombination  probabilities  calculated  for 
different  adsorption  geometries  of  the  two  ad- atoms  (see  Ref.  [21]  for  details).  The  L-H  simulation  was 
performed  at  three  different  temperatures  of  the  substrate:  Ts=600K,800K,1000K. 

The  semiclassical  recombination  coefficient  y  for  the  E-R  and  the  L-H  recombination  is  reported  in  Fig. 
10,  together  with  the  high  temperature  experimental  data  determined  recently  by  M.Balat  et  coll,  on  (3- 
cristobalite  using  the  MESOX  set  up  [31,32]. 
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Fig.  10  :  The  semiclassical  recombination  coefficient  y  for  the  E-R  (X )  and  the  L-H  (+) 
recombination  reaction  [21]  is  reported  together  with  the  experimental  data 
(  aScott  (silica-based  materials)Ref.  [3];  •  Greaves  et  al.  (silica)  Ref.  [4]  ;  Berkut  et  al. 
(silica-based  materials)  Ref.  [2]  ;  0  Balat  et  al.  (p-cristobalite)  Ref.  [31] ). 


A  cross-comparison  between  the  experimental  and  the  semiclassical  results  shows  that  a  good  agreement 
is  reached  between  the  E-R  y  values  and  the  experimental  values,  whereas  at  the  lower  temperature  the  L- 
H  recombination  probabilities  match  the  experimental  measurements.  Note  that  a  large  surface 
temperature  effect  is  found  for  the  L-H  recombination  process,  whereas  Ts  has  only  a  weak  impact  on  the 
E-R  reaction. 

The  theoretical  y  values  should  not  be  directly  compared  with  the  cell  measurements  (see  the  previous 
discussion  on  this  point),  nevertheless  the  agreement  between  theory  and  experiments  is  very  encouraging. 

The  reliability  of  the  semiclassical  approach  followed  in  the  MD  simulation  has  been  demonstrated  in  a 
combined  experimental  and  theoretical  investigation  [22]  recently  carried  out  on  the  catalytic  activity  of 
the  P-quartz.  In  the  table  below  the  theoretical  and  experimental  recombination  coefficient  y  for  oxygen 
recombination  on  P-quartz  are  reported  at  TS=1000K  and  compared  with  the  corresponding  data  for  P- 
cristobalite. 
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Surface  (Ts=1000K) 

Yth 

Yexp 

P-quartz 

0.008 

0.008 

P-cristobalite 

0.03 

0.027 

Such  agreement  gives  support  to  the  E-R  recombination  mechanism  as  the  mechanism  followed  in  the 
recombination  reaction  process  monitored  in  the  experiments  and  gives  further  evidence  of  the  large  effect 
the  surface  crystallographic  structure  has  on  the  catalytic  activity  of  silica-based  surfaces. 

It  is  worth  noticing  that  the  semiclassical  recombination  coefficients  reported  in  the  table  above  were 
obtained  assuming  that  in  the  initial  conditions  the  O  ad-atom  is  randomly  adsorbed  within  the  unit  cell. 
The  same  trend  is  observed  when  the  O  ad-atom  is  chemisorbed  in  a  specific  silica  site,  on  top  of  a  Si 
surface  atom.  In  this  case  in  fact  we  get  y=  0.04  and  y=0.3  for  O  recombination  on  P-quartz  and  P- 
cristobalite  at  Ts=1000K,  respectively.  That  is,  although  the  recombination  coefficients  are  higher  than 
that  measured  in  the  experiments,  the  ratio  y(cristobalite)/y(quartz)  is  about  a  factor  two  higher  compared 
to  the  ratio  observed  in  the  experiments. 

The  semiclassical  MD  calculations  have  revealed  some  aspects  that  cannot  be  observed  directly  in  the 
experiments.  The  dynamics  is,  in  fact,  quite  complex  with  several  collisional  processes  that  are  in 
competition  with  the  O2  formation  at  the  surface. 


O  +  Oads  *silica 

- >  02(v,  j)+silica 

recombination 

(21) 

Ogas  "T  Ogas  +silica 

inelastic 

Oads  *silica+  Ogas 

Ogas  "F  Oads  *  silica 

adsorption/desorption 

[02]ads 

Oads  +  Oads  +  silica 

adsorption 

In  particular,  the  adsorption  and  adsorption/desorption  processes  of  the  oxygen  atoms  are  quite  effective  in 
a  large  collisional  energy  range.  This  is  shown  in  Fig.  11  where  the  probabilities  for  the  different  surface 
processes  active  on  P-cristobalite  are  reported  at  TS=1000K. 
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Kinetic  Energy(eV) 


Fig.  11 :  Probability  for  the  different  surface  processes  of  Eq.  (21)  on  p-cristobalite.  Ts=1000K. 

In  the  full  range  of  the  explored  energies  the  adsorption  and  adsorption/desorption  processes  are  quite 
effective.  This  is  an  interesting  result  that  could  be  of  some  relevance  for  the  surface  damage  problem.  The 
residence  time  for  these  processes  is  large  so  that  chemical  reactions  between  the  adsorbed  atoms  and  the 
surface  Si/O  atoms  can  be  eventually  be  promoted  (with  the  formation  of  SiOx  volatile  compounds). 

In  Fig.  12  the  probability  for  the  different  reaction  channels  of  Eq.  (21)  is  reported  but  for  the  interaction 
of  N  atoms  with  P-cristobalite.  In  the  initial  conditions  the  N  ad-atom  is  chemisorbed  in  a  specific 
chemisorption  site,  that  is  on  top  of  a  Si  surface  atom,  in  thermal  equilibrium  with  the  silica  surface  at 
Ts=1000K.  The  gas-phase  N  hits  the  surface  in  the  normal  configuration,  with  the  (X,Y)  coordinates 
randomly  varied  within  the  unit  silica  cell  of  the  P-cristobalite.  The  potential  energy  surface  assumed  in 
the  dynamics  is  the  LEPS  surface  of  Fig.  4  obtained  from  DFT  calculations.  A  further  aspect  of  great 
relevance  for  the  surface  heating  problem  concerns  the  energetics  of  the  recombination  process.  Reactions 
(18),  (19)  are  both  strongly  exothermic  so  that  part  of  the  exothermic  energy  can  be  gained  by  the  formed 
02  molecules  while  the  remaining  part  can  be  transferred  to  the  substrate  as  heat  flux.  The  energy  flux 
pathways  are  shown  in  Fig  13.  Here  the  energy  partitioning  among  the  translational  (Etr),  rotational  (Erot) 
and  vibrational  motions  (Evib)  of  the  02  molecules  and  the  surface  phonons  (AEph)  is  shown  as  a  function 
of  the  impact  kinetic  energy  for  oxygen  formation  on  P-cristobalite  at  Ts=1000K  (a)  and  Ts=1600K, 
respectively. 
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Fig.  12:  The  same  as  for  Fig. 11  but  for  N  +  Nad  on  (3-cristobalite.  Nad  is  initially  adsorbed 
on  top  of  a  Si  surface  atom.  The  PES  of  Fig.  4  was  assumed  in  the  dynamics. 


13a 


13b 


Fig.:  13:  Energy  partitioning  among  the  translational  (Etr),  rotational  (Erot),  vibrational  (EVib) 
motions  of  the  O2  molecules  and  the  surface  phonons  (AEph)  as  a  function  of  the  impact  kinetic 
energy  in  the  reaction:  O  +  Oad  *silica  --->  02(v,j)  +  silica  .  (a)  Ts=1000K  and  (b)  Ts=1600K. 


It  is  evident  that  the  largest  fraction  is  primarily  shared  among  the  vibrational  and  rotational  motion  of  02, 
while  only  a  small  fraction  of  the  total  energy  is  transferred  to  the  silica  surface  as  phonon  excitation 
energy.  However,  at  the  higher  energies  and  Ts=1000K  almost  30%  of  the  exothermic  energy  appears  as 
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phonon  energy.  Interesting  different  behaviours  of  the  reaction  energetics  can  be  pointed  out  by  comparing 
the  results  reported  at  the  two  different  surface  temperatures.  As  a  consequence  of  the  energy  sharing 
mechanism,  the  02  molecules  are  formed  in  highly  vibrationally  excited  state. 

This  is  shown  in  Fig.  14  where  a  typical  vibrational  distribution  is  reported  for  the  recombination  reaction 
over  (3-cristobalite  at  Ekin=  0.003eV,  Ts=1000K. 


Fig. 14:  A  typical  vibrational  distribution  of  the  O2  molecules  formed  after 
atom  recombination  over  p-cristobalite  at  Ekm=  0.003eV  and  TS=1000K. 


The  distributions  are  characterized  by  a  peak  in  v=0,  but  virtually  all  the  vibrational  levels  are  populated 
and  involved  in  the  surface  process.  The  distributions  are  in  non-thermal  equilibrium. 

From  the  calculated  vibrational  distribution,  the  vibrational  state- selected  recombination  coefficient 
y(v,Ts)  can  be  easily  obtained.  The  spectral  distribution  of  y  for  O  recombining  on  P-cristobalite  at 
Ts=1000K  is  shown  in  the  Fig.  15.  The  y(v)  coefficients  are  an  important  ingredient  for  the  kinetic 
modelling  of  the  chemistry  of  the  boundary  layer  [33-34]  where,  very  likely,  a  non-thermal  equilibrium 
condition  is  established  due  to  the  fluodynamics  conditions. 
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Fig. 15:  The  vibrational  spectral  distribution  of  the  recombination  coefficient  y  for  02(v) 
formed  in  the  reaction:  O  +  Oad  *silica  --->  02(v,j)  +  silica  is  reported  as 
a  function  of  the  vibrational  quantum  number  v.  Ts=1000K. 

The  energy  transferred  to  the  surface  in  the  recombination  processes  is  generally  claimed  to  be  among  the 
primary  sources  of  the  surface  damage.  Molecular  dynamics  simulations  for  O  (  and  N)  recombination  on 
silica  surfaces  have  shown  that  this  energy  is  only  a  small  fraction  of  the  total  energy  delivered  in  the 
reaction,  so  that  the  energy  accommodation  coefficient  P  should  be  less  than  one  (  contrary  to  what  is 
generally  assumed  in  kinetic  modelling  calculations,  i.e  P=1  ).  Indeed,  the  overall  energy  transferred  to  the 
substrate  is  the  result  of  a  complex  mechanism  and  it  is  due  to  the  concurrence  of  various  surface 
processes.  The  adsorption  processes  can  play  an  important  role  since  the  energy  exchanged  with  the 
surface  can  be  quite  large.  In  Fig.  16  a  typical  adsorption  collision: 

O  +  Oad  +  silica - >  [02]ad  +  silica - >  O  +  Oad  +  silica/(3-cristobalite  Ts=1000K 

is  reported  together  with  the  corresponding  energy  transferred  to  the  silica  substrate. 
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Fig.  16:  A  typical  collision  for  the  adsorption  process:  O  +  Oad  +  silica  — >[02]ad  +  silica— > 
O  +  Oad  +  silica/p-cristobalite  is  reported  together  with  the  corresponding  energy 
transferred  to  the  silica  substrate,  a)  the  z  coordinate  of  the  two  oxygen  atoms 
is  reported  as  a  function  of  the  collisional  time,  b)  the  corresponding 
energy  transferred  to  the  silica  surface  is  shown.  Ts=1000K. 


FINAL  REMARKS 

The  development  and  applications  of  the  semiclassical  approach  to  the  molecule- surface  dynamics  have 
made  possible  MD  simulations  of  elementary  surface  processes  involving  oxygen  and  nitrogen  atoms  over 
silica  model  surfaces  relevant  to  the  TPS  system.  In  the  simulations  various  important  behaviours  of  the 
surface  substrate,  such  as  the  phonon  structure,  the  surface  polymorphic  modifications,  the  site  and  the 
intrinsic  surface  temperature  effect  can  be  incorporate  and  their  impact  on  the  surface  catalytic  activity 
studied  at  a  very  elementary  level.  Collisional  data,  such  as  the  global  and  the  vibrational  state  selected 
recombination  coefficient  y,  the  energy  distributions  in  the  product  states,  the  energy  flow  pathways  and 
the  energy  released  to  the  surface  in  the  recombination  reaction  as  well  as  in  the  adsorption/desorption 
processes  can  be  determined.  The  accuracy  of  the  calculated  coefficients  mostly  relies  on  the  accuracy  of 
the  interaction  potential  assumed  in  the  dynamics,  rather  than  to  the  approximations  made  in  the 
collisional  model.  As  a  matter  of  fact,  the  determination  of  a  reasonable  PES  for  these  complex 
heterogeneous  systems  is  the  most  critical  problem  to  solve.  MD  simulations  have  some  limitations, 
among  the  others:  i)  the  surface  processes  can  be  followed  in  the  short  time  scale  (few  picoseconds).  ii)the 
surface  inserted  in  the  simulations  is,  in  general,  a  flat  perfect  crystal  surface,  iii)  the  dynamics  occurs  on 
an  adiabatic  potential  surface  so  that  surface  processes  involving  electronically  excited  states  are  not 
explored.  Though  the  inclusion  of  these  effects  add  more  complexity  and  computational  time 
requirements,  the  importance  of  MD  calculations  for  the  catalytic  systems  candidate  for  the  future  TPS 
systems  will  become  more  evident  in  the  next  future. 


5-20 


RTO-EN-AVT-1 42 


1 u 

NATO 

3! 

OTAN 

Molecular  Dynamics  Simulations  of  Surface  Processes: 
Oxygen  Recombination  on  Silica  Surfaces  at  High  Temperature 


REFERENCES 

[1]  M.C.Lin  and  GJ.Ertl  Ann.Rev.Phys.Chem.  37  (1986)  587;  G.J.Ertl  Catal.Rev.-Sci.Eng.  21  (1980) 
201 

[2]  V.D.  Berkut,  V.M.  Doroshenko,  V.V.  Kovtun,  N.N.  Koudryavtsev,  S.S.  Novikov,  N.V.  Smirnov, 
A.I.  Sharotovov,  Soviet  J.  Chem.  Phys.  9  (1992)  2222 

[3]  C.D.  Scott,  AIAA  Paper  83-1485  (1983). 

[4]  J.C.  Greaves,  J.W.  Linnett,  Trans.  Faraday  Soc.,  55  (1959)  1355 

[5]  E.  J.  Jumper,  W.A.  Seward,  J.  Thermophys.  &  Heat  Transfer  8  (1994  )  460 

[6]  F.  Nasuti,  M.  Barbato,  C.  Bruno,  J.  Thermophysics  &  Heat  Transfer  10  (1996)  131 

[7]  V.  L.  Kovalev,  “ Phenomenological  treatment  of  heterogeneous  recombination”,  these  VKI  Lecture 
Notes 

[8]  A.E.Depristo  and  A.Kara  Adv.Chem.Phys.  77  (1990)  163 

[9]  G.P.Brivio  and  T.B.Grimley  Surf.Sci.Rep. 17(1993)  1 

[10]  A.Gro() :  these  VKI  Lecture  Notes 

[11]  R.  Dovesi,  et  al  CRYSTAL  95  User's  Manual’,  (  University  of  Torino,  Torino,  1996). 

[12]  E.g.,  see,  K.  G.  Nakamura,  Chem.  Phys.  Lett.  285  (1998)  21. 

[13]  W.  Koch,  W.G.  Holthousen,  ‘A  chemist’s  guide  to  density  functional  theory’,  (Wiley-VCH, 
Weinheim,  2000). 

[14]  14)  M.  Rutigliano,  A.  Pieretti,  M.  Cacciatore,  N.  Sanna  and  V.  Barone,  “N  atoms  recombination  on 
a  silica  surface:  a  global  theoretical  approach”.  Surf.  Sci.  2006 

[15]  C.  Lee,  W.  Yag  and  R.  G.  Parr,  Phys.  Rev.  37  (1988)  785. 

[16]  Gaussian  03  (Revision  C.02),  M.  J.  Frisch,  et  al,  Gaussian,  Inc.,  Wallingford  CT,  2004. 

[17]  “Spectra  of  Diatomic  Molecules”,  Van  Nostrand,  Princeton,  NJ,  1950;  K.P.  Huber  and  G.  Herzberg, 
“  Molecular  Spectra  and  Molecular  Structure,  IV.  Constants  of  Diatomic  Molecules”,  Van  Nostrand, 
New  York,  1979. 

[18]  J.  McCreery  and  G.  Walken,  J.  Chem.  Phys.  63  (1975)  2340 

[19]  B.  Jackson  and  M.  Persson,  J.  Chem.  Phys.  96  (1992)  2378 

[20]  M.Cacciatore,  E.Christoffersen,  M.Rutigliano  J.  Phys.  Chem.  A,  108  (2004)  8810 

[21]  M.Cacciatore, M.Rutigliano  and  G.D.Billing  ,  J.  Therm, &Heat  Transfer.  13  (1999)195 

[22]  L.  Bedra,  M.  Rutigliano,  M.  Balat-Pichelin,  M.  Cacciatore,  “Atomic  Oxygen  Recombination  on 
Quartz  at  High  Temperature:  Experiments  and  MD  Simulation”,  submitted 


RTO-EN-AVT-1 42 


5-21 


Molecular  Dynamics  Simulations  of  Surface  Processes: 

Oxygen  Recombination  on  Silica  Surfaces  at  High  Temperature 


[23]  Van  Beest  B.  W.  H„  Kramer  G.  J.,  Van  Santen  R.  A.,  Phys.  Rev.Lett.,  64  (1990)1955. 

[24]  G.D.Billing  "Dynamics  of  Molecule  Surface  Interactions",  (Wiley  &  Sons,  N.  Y.,  2000). 

[25]  G.  D.  Billing  and  M.  Cacciatore,  Chem.  Phys.  Lett.  113  (1985)  23 

[26]  G.D.Billing,  M.Cacciatore  Chem.Phys.103  (1986)  137 

[27]  S.Iannotta,  C.Gravili,A.Boschetti,A.Cagol  and  M.Cacciatore  Chem.Phys.  194  (1995)133 

[28]  R.W.G.Wickoff  American  Journal  Science  448  (1925)  5 

[29]  B.P.  Feuston  and  S.  H.  Garofalini,  J.  Chem.  Phys.  91  (1989)  564 

[30]  M.  Balat-Pichelin,  J.  M.  Badie,  M.Cacciatore,  M.Rutigliano,  ESA  SP-521  (2003)  403 

[31]  M.  Balat-Pichelin,  J.M.  Badie,  R.  Berjoan,  P.  Boubert,  Chem.  Phys.  291  (2003)  181 

[32]  M.Balat-Pichelin  these  VKI  Lecture  Notes 

[33]  I.  Armenise,  M.  Capitelli,  C.  Gorse,  M.  Cacciatore  and  M.  Rutigliano,  JSR  37  (2000)  318 

[34]  D.  Bruno,  M.  Cacciatore,  S.  Longo,  M.  Rutigliano,  Chem.  Phys.  Lett.  320  (2000)  245 


5-22 


RTO-EN-AVT-142 


